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We review several topics related to the diagonalization of quantum field Hamilto- 
nians using the quasi-sparse eigenvector (QSE) method. 



1 Introduction 

Quasi-sparse eigenvector (QSE) diagonalization is a new computational 
method which finds the low-lying eigenvalues and eigenvectors for a general 
quantum field Hamiltonian . It is able to handle the exponential increase 
in the size of Fock space for large systems by exploiting the sparsity of the 
Hamiltonian. QSE diagonalization can even be applied directly to infinite- 
dimensional systems. The method is most effective when the splitting be- 
tween low-lying eigenvalues is not too small compared to the size of the off- 
diagonal Hamiltonian matrix entries. In such cases the low- lying eigenvectors 
are quasi-sparse, meaning that the vector is dominated by a small fraction of 
its largest components. The QSE algorithm can then be applied to the 
Hamiltonian H using the following steps: 

1. Select a subset of orthonormal basis vectors {e^, • • • ,e, n } and call the 
corresponding subspace S. 

2. Diagonalize H restricted to S and find one eigenvector v. 

3. Sort the basis components (e^ \v) according to their magnitude and re- 
move the least important basis vectors. 

4. Replace the discarded basis vectors by new basis vectors. These are 
selected at random from a pool of candidate basis vectors which are 
connected to the old basis vectors through non- vanishing matrix elements 
of if. 

5. Redefine S as the subspace spanned by the updated set of basis vectors 
and repeat steps 2 through 5. 

If the subset of basis vectors is sufficiently large, the exact eigenvectors will 
be stable fixed points of the update process. 
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The purpose of this short tutorial is to provide additional information 

for those interested in writing or using codes that implement QSE methods. 

We discuss the calculation of Hamiltonian matrix elements in the momentum 

Fock space representation and the generation and selection of new basis states. 

We also include a simple program which applies the QSE method to find the 

ground state of an input matrix. Readers interested in an introductory 

overview of QSE diagonalization are encouraged to consult the references in 
1 



2 Fock states and the Hamiltonian matrix 

It is most effective to describe the details of a QSE calculation in the context 
of a specific example. In this discussion we consider <fi 4 theory in 1 + 1 
dimensions, the same example as in . Some of the techniques described 
here have been specifically optimized for scalar field theories. Other systems 
require a somewhat different set of tools and will be discussed in the future. 
The Hamiltonian density for </> 4 theory in 1 + 1 dimensions has the form 

" 2 + A. 



rL ~ 2 \dt J ~r 2 \dx J ^2 



2 V ~ IT-V -J 

where the normal ordering is with respect to the mass /U. We consider the 
system in a periodic box of length 2L. We then expand in momentum modes 
and reinterpret the problem as an equivalent Schrodingcr equation 2 . We 
will use a momentum cutoff scheme where the modes q n are restricted in 
momentum so that \n\ < N max . The resulting Hamiltonian is 

H = -2-EWbal: + IE faM - «-»*• (!) 

n n 
ni +n 2 +n s +n 4 =Q 

where 



^)^TT+I' 2 (2) 
and b(fi) is the coefficient for the mass counterterm 

n 

It is convenient to split the Hamiltonian into free and interacting parts 
with respect to an arbitrary mass fjl: 
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n 



(4) 
(5) 



<Zni <Zn 2 9"3 9"4 ' 



11+12+13+14=0 



// will be used to define the basis states of our Fock space. Since H is 
independent \j! , it is useful to perform calculations for different /j,' to obtain 
an estimate of the error. It is also helpful to find the range of values for // 
which maximizes the quasi-sparsity of the eigenvectors. This can be used to 
improve the accuracy of the QSE calculation. 

Let us define ladder operators with respect to 



a„(ju') 



1 



These satisfy the usual commutation relations, 

[a n ,a m ] = [al,aln] = °- 
The Hamiltonian can now be rewritten as 



(6) 
(7) 



(8) 
(9) 



(10) 



+ 192L 

"1 +"2 +"3+14 = 



( a "l+ a -» 1 ) ( a "2+ a -n 2 ) ( a "3+ a - re3 ) ( a "4+ a -„ 4 ) 
\f^JP) \f^U^) y/^^ifJ-') y/ul^ifJ,') 



In (10) we have omitted constants contributing only to the vacuum energy. 

We can represent any momentum-space Fock state as a string of occupa- 
tion numbers, 



°-N„ 



.OAT max ) 



(11) 



where 



ala n \o- Nrr 



.°JV m «> = O n |0_ 



1 OiVmax) • 



(12) 
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From the usual ladder operator relations, it is straightforward to calculate 
the matrix element of H between two arbitrary Fock states. In 1 auxiliary 
cutoffs were used to render the dimension of Fock space finite. This is in fact 
an unnecessary restriction and QSE diagonalization encounters no difficulties 
dealing with the full infinite-dimensional space. 

3 New basis vectors 

Aside from calculating matrix elements, the only other fundamental operation 
involving basis vectors is the generation of new basis vectors. As mentioned 
before, the new states should be connected to an old basis vector through non- 
vanishing matrix elements of H . Let us refer to the old basis vector as |e). 
For this example there are two types of terms in our interaction Hamiltonian, 
a quartic interaction 

T3 ( a "i ~l~ a -7ii) ( a «2 + a -n 2 ) ( a "3 a -™;?) (^ a -ni-n 2 —n 3 + a ni+n 2 +n 3 ) ; 

ni, 712,^3 

(13) 

and a quadratic interaction 

53 ( a -« + 4) (<»n + at n ) • ( 14 ) 

n 

To produce a new vector from |e) we simply choose one of the possible operator 
monomials 

Q j ni®ri2®n3Q'—ni—n2—n3i <n\ ®n2 ^713^— n\ — n 2 — n$ > ' ' > 

and let it act upon \e). Our experience is that the interactions involving 
the small momentum modes are generally more important than those for the 
large momentum modes, a signal that the ultraviolet divergences have been 
properly renormalized. For this reason it is best to arrange the selection 
probabilities such that the smaller values of |m|, \n 2 \, \n 3 \ and \n\ are chosen 
more often. 

We note that the new basis vector selection will occasionally fail. Either 
the vector is zero due to an annihilation operator acting on an unoccupied 
state or the vector is already in our subset of basis vectors and therefore not 
new. In either case we simply select a new basis vector again. If the selection 
process fails repeatedly then a different old basis vector |e') is used. 
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(15) 



4 Sample code 



We have considered two basic operations, the calculation of matrix elements 
and the generation of new basis states. The subroutines which perform these 
tasks will depend on the form of the quantum held Hamiltonian. The main 
structure of the QSE algorithm, however, is independent of the details of the 
Hamiltonian. We demonstrate the essential features of the algorithm with 
the following short MATLAB program. In this example we find the ground 
state of a finite matrix. To keep the example as simple as possible, we avoid 
calculating matrix elements and generating new basis states by assuming that 
the entire Hamiltonian matrix can be stored in memory. 

niter = 30; 
nactive = 100; 
nretain = 80; 

H = loadsparse ( ' samplematrix' ) ; 
vecs = [1 inactive]; 
for iter = l:niter 

[v,d] = eigs (H(vecs ,vecs) , 'SR' ,1) ; 

weights!. = v . ~2 ; 

[sortelements , sortorder] = sort (-weightsl) ; 
vecs = vecs (sortorder (1 :nretain) ) ; 
wtsuml = cumsum (weightsl) /sum(weightsl) ; 
while (length(vecs) <nactive) 

overl = f ind(rand*wtsuml (nretain) <wtsuml) ; 

seedvec = vecs (overl (1) ) ; 

candidates = find(H( :, seedvec) ) ; 

candidates = candidates (candidates ~= seedvec); 

weights2 = abs (H(candidates , seedvec) ) ; 

wtsum2 = cumsum (weights2) /sum(weights2) ; 

over2 = find(rand <wtsum2) ; 

newvec = candidates (over2 (1) ) ; 

vecs = unique ([vecs newvec]); 

end 

end 

save results. mat vecs v d 

niter is the number of iterations of the algorithm, nactive is the number 
of basis states used to form the subspace S, and nretain is the number 
of basis states kept after removing the tail of the eigenvector distribution, 
loadsparse is a program which reads in a datafile called samplematrix. This 
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data is stored in the matrix variable H. The remaining code can be broken 
down in terms of the five basic steps of the QSE algorithm. 

Select a subset of orthonormal basis vectors {e il7 ■ ■ ■ , e in } and call the corre- 
sponding subspace S. 

vecs = [1 inactive] ; 
The variable vecs is a row vector which labels the basis vectors in the subspace 
S. In our example vecs is initialized to be the first nactive basis vectors. We 
are assuming that the ground state of the system has a non-zero overlap with 
at least one of the first nactive basis vectors. A more general initialization 
would be to choose a random subset of basis vectors. 

Diagonalize H restricted to S and find one eigenvector v. 

[v,d] = eigs (H(vecs ,vecs) , 'SR' ,1) ; 
The function eigs is a sparse matrix Arnoldi diagonalization routine. It is 
being called with the parameters 'SR' ,1 which tells eigs to find the eigen- 
value and eigenvector with the smallest real part. The Hamiltonian submatrix 
corresponding with basis vectors in vecs is diagonalized. The eigenvector is 
stored in the column vector v and the eigenvalue is stored as the variable d. 

Sort the basis components (e ij \v) according to their magnitude and remove 
the least important basis vectors. 
weightsl = v. ~2; 

[sortelements, sortorder] = sort (-weightsl) ; 

vecs = vecs (sortorder (1 :nretain) ) ; 
weightsl is a column vector which contains the squares of the components 
of the eigenvector v. sortelements is a list of the elements of -weightsl 
from smallest to greatest, sortorder is the corresponding list of locations 
for these elements. If for example weightsl = [ . 2 .5 . 3] , then sortorder 
= [2 3 1]. The list sortelements is not used. Only the nretain most 
important basis vectors are kept in vecs. 

Replace the discarded basis vectors by new basis vectors. These are selected 
at random from a pool of candidate basis vectors which are connected to the 
old basis vectors through non-vanishing matrix elements of H. 
wtsuml = cumsum(weightsl) /sum(weightsl) ; 
while (length(vecs) <nactive) 

overl = find(rand*wtsuml (nretain) <wtsuml) ; 

seedvec = vecs (overl (1) ) ; 

candidates = f ind(H( :, seedvec) ) ; 

candidates = candidates (candidates ~= seedvec); 

weights2 = abs (H(candidates , seedvec) ) ; 
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wtsum2 = cumsum(weights2) /sum(weights2) ; 
over2 = f ind(rand <wtsum2) ; 
newvec = candidates (over2 (1) ) ; 
vecs = unique ( [vecs newvec]); 

end 

cumsum is a function which takes the input list and creates another list con- 
sisting of partial sums, wtsuml is therefore an ascending string of numbers 
between and 1 such that the difference between an entry and its preceding 
entry is proportional to the corresponding value of weights 1. overl is a 
location list of elements of wtsuml which are greater than the product of a 
generated random number between and 1 and wtsuml (nretain) . overl (1) 
is the first entry of overl. We note that overl (1) must be less than or equal 
to nretain. seedvec is the basis vector which corresponds with overl (1). 
candidates is a column vector consisting of all basis vectors which have a 
non-zero matrix element with seedvec. Since the task is to find new basis 
vectors, seedvec is removed from candidates. 

weights2 is a column vector containing the absolute values of the matrix 
transition elements from seedvec. wtsum2 is an ascending string between 
and 1 such that the difference between an entry and its preceding entry is 
proportional to the value of weights2. over2 is a location list of elements of 
wtsum2 greater than a generated random number, newvec is the basis vector 
which corresponds with over2(l) . newvec is appended to the set vecs. The 
function unique sorts the input list and removes any repetitions. This loop 
is repeated until vecs has nactive elements. 

Redefine S as the subspace spanned by the updated set of basis vectors and 
repeat steps 2 through 5. 
for iter = l:niter 

end 

save results. mat vecs v d 
The algorithm is iterated niter times and the final results are saved in the 
file results. mat. More examples of QSE diagonalization codes for various 
field theory systems written in MATLAB, Fortran, and/or C++ will be made 
available in the future. 

5 Summary 

In this short tutorial we have listed some supplemental material for those 
interested in writing or using codes that implement QSE methods. We have 
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discussed the calculation of Hamiltonian matrix elements and the selection of 
new basis states. We have also presented a simple MATLAB program which 
applies the QSE method to find the ground state of a matrix. 

The discussion here was limited to the basic QSE method. As presented 
in the Guangzhou workshop there have also been interesting new developments 
in QSE diagonalization regarding the technique of stochastic error correction. 
This technique will be presented in a forthcoming work 4 . 
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